clear; clc; close all;
% load XYSVLambda2_LD
suffixes = ['0', '1', '2', '3', '4'];
betaCols = [1,2,3,4,5,2758,2759,2760,2761,2762,2763,2764,2765,2766,2767,2768,2769,2770,2771,2772,2773,2774,125878,125879,125880,125881,125882,125883,125884,125885,125886,125887,125888,125889,125890,125891,125892,125893,125894,125895,125896,125897,125898,125899,125900,125901,125902,125903,125904,125905,125906,125907,125908,125909,125910,125911,125912,125913,125914,125915,125916,125917,125918,];
for suffix = suffixes,
    clearvars -except suffix suffixes betaCols
    data_MM = load(['C:\Users\Mittal\Documents\Work\Code\survival-analysis\perl\NTDBNEW_Imputation_Train_LD_imputed_' suffix '_MM.txt']);
    data = spconvert(data_MM);
    data = data(:,[1 (betaCols+1)]);
    Y = data(:,1);
    X = data(:,2:end);
    
    thresh = 1e-7;
    
    fid = fopen(['NormalizedX_LD_' suffix '.txt'],'w');
    count = 1;
    for i = 1:size(X,2),
        numNonzero(i) = nnz(X(:,i));
        D = X(:,i);
        fprintf('%d i = %d\n',count,i);
        if(nnz(X(:,i)) > 1534)
            D = D - mean(D);
            count = count + 1;
        end
        if(norm(D) ~= 0)
            D = D/norm(D);
            numSmall = full(sum(abs(D) < thresh));
            numZero = full(sum(abs(D) == 0));
            if(numSmall-numZero > 0)
                fprintf('Setting %d entries to 0\n', numSmall-numZero);
                D(abs(D) < thresh) = 0;
            end
            if(mod(i,1000) == 0)
                fprintf('%d i = %d\n',count,i);
            end
            [j,k,val] = find(D);
            for ii = 1:length(j),
                fprintf(fid,'%d %d %d\n',j(ii),i,val(ii));
            end
        end
    end
    fclose(fid);
    clear X;
    
    X = load(['NormalizedX_LD_' suffix '.txt']);
    X = spconvert(X);
    
    [~,S,V] = svds(X);
    
    fid = fopen(['Lambda2_LD_' suffix '.txt'],'w');
    for i = 1:size(V,1),
        %     if(mod(i,1000) == 0)
        fprintf('i = %d\n',i);
        %     end
        lambda2Temp = sparse((V(i,:).*diag(S)'*V').^2);
        lambda2Temp(abs(lambda2Temp) < thresh) = 0;
        [j,k,val] = find(lambda2Temp);
        for ii = 1:length(j),
            fprintf(fid,'%d %d %d\n',i,k(ii),val(ii));
        end
    end
    fclose(fid);
    
    lambda2Temp = load(['Lambda2_LD_' suffix '.txt']);
    lambda2 = spconvert(lambda2Temp);
    if(lambda2Temp(end,1) < size(V,1))
        lambda2 = [lambda2; zeros((size(V,1)-lambda2Temp(end,1)),lambda2Temp(end,1))];
        lambda2 = [lambda2 zeros(size(V,1),1)];
    end

    save(['XYSVLambda2_LD_' suffix], 'X', 'Y', 'S', 'V', 'thresh', 'lambda2', 'data');
    betaTemp = ReadModelFile(['modelfile_LD_L1_' suffix  '.txt']);
    intercept = betaTemp(1,2);
    betaTemp = sortrows(betaTemp(2:end,:));
    beta = zeros(length(betaCols),1);
    for i = 1:length(beta),
        if(find(betaTemp(:,1) == betaCols(i)))
            beta(i) = betaTemp(find(betaTemp(:,1) == betaCols(i)),2);
        end
    end
    
    betaX = intercept + data(:,2:end)*beta;
    LPredY = 1./(1 + exp(-betaX));
    lYhat = log(LPredY./(1-LPredY));
    stdlYhat = std(lYhat);
    Rsq = (sum((LPredY-mean(Y)).^2))/(sum((Y-mean(Y)).^2));
    betaStand = beta*(sqrt(Rsq)/stdlYhat);
    epsilon = lambda2*(betaStand.^2);
    PropWeights = epsilon/sum(epsilon);
    
    RawRelWeight = epsilon;
    RescaledRelWeight = PropWeights;
    keys = textread('NTDBNEW_Imputation_Keys.txt','%s','whitespace','\n');
    keys = keys(betaCols);
    fid = fopen(['RawRelWeights_LD_L1_' suffix '.txt'],'w');
    for i = 1:length(RawRelWeight),
        fprintf(fid,'%s#$#$#%d\n',keys{i}, RawRelWeight(i));
    end
    fclose(fid);
    
    fid = fopen(['RescaledRelWeights_LD_L1_' suffix '.txt'],'w');
    for i = 1:length(RescaledRelWeight),
        fprintf(fid,'%s#$#$#%d\n',keys{i}, RescaledRelWeight(i));
    end
    fclose(fid);
    
    fid = fopen(['beta_LD_L1_' suffix '.txt'],'w');
    for i = 1:length(beta),
        fprintf(fid,'%s#$#$#%d\n',keys{i}, beta(i));
    end
    fclose(fid);
end